# Back of the envelope calculations of damages 
library(pacman)
p_load(
  here, fst, data.table, ggplot2, magrittr, janitor, stringr, purrr
)

# Loading data ----------------------------------------------------------------
  # Marginal production model
  model_dt = 
    #read.fst(
    #  path = here("Data/electricity-generation/plant-model-fit/good-model-dt.fst"),
    #  path = here("Data/electricity-generation/plant-model-fit/own-region-model-dt.fst"),
    #  as.data.table = TRUE
    #)[,plant_id_eia := as.character(plant_id_eia)] 
   fread(
     here("Data/electricity-generation/output/PPRegressors-oge.csv")
    )[,plant_id_eia := as.character(orispl_code)] |> 
    setkey(plant_id_eia)
  setnames(
    model_dt, 
    old = str_subset(colnames(model_dt), 'sq'),
    new = paste0(
      'excess_load_sq_',
      str_subset(colnames(model_dt), 'sq') |>
        str_remove('_sq') |>
        str_remove('excess_load_')
  ))
  setnames(
    model_dt, 
    old = str_subset(colnames(model_dt), 'excess_load'),
    new = paste0('coef_',str_subset(colnames(model_dt), 'excess_load'))
  )
  # Load data
  nerc_load_dt = 
    fread(
      here("Data/electricity-generation/output/RegLoad-oge.csv")
    )[,.(
      nerc_adj, 
      datetime_utc = utc_time, 
      excess_load, 
      excess_load_sq = excess_load^2
    )] |>
    dcast(
      formula = datetime_utc ~ nerc_adj,
      value.var = c('excess_load','excess_load_sq')
    ) |> 
    clean_names() |>
    setkey(datetime_utc) 
  # Total Load
  
  # Marginal damages 
  tot_md_dt = fread(
    here('Data/electricity-generation/output/damage-per-mwh-oge.csv')
  )[,-'V1']
  # Avg solar panel production by region (From Mark email 10/17/2022)
  panel_output_dt = data.table(
    nerc_adj = c('cal','mro','npcc','rfc','serc','tre','wecc'),	
    avg_output_mwh = c(5935.968, 4802.954, 4546.329, 4584.62, 5161.415, 5274.348, 5444.394)/1000
  )
  # Time profile for solar production by state
  state_solar_profile_dt = fread(
    here("Data/electricity-generation/output/solar-profile-dt.csv")
  )
  # Crosswalk between states and nerc regions 
  nerc_state_xwalk = data.table(
    nerc_adj = c('cal','mro','npcc','rfc','serc','tre','wecc'),
    state_fips = c(6, 29, 25, 39, 13, 48, 49)
  )
  # Plant info table
  oge_plant_info_dt =
    rbind(
      read.fst(
        path = here("Data/electricity-generation/plant-info-dt-oge.fst"),
        as.data.table = TRUE
      )[, 
        .SD[which.min(year)], 
        by = .(plant_id_eia)
      ],
      read.fst(
        path = here("Data/electricity-generation/shaped-info-dt-oge.fst"),
        as.data.table = TRUE
      )[, 
        .SD[which.min(year)], 
        by = .(plant_id_eia)
      ],
      fill = TRUE
    )[plant_id_eia %in% model_dt$plant_id_eia,
      .(plant_id_eia = as.character(plant_id_eia), namepcap, stack_height)
    ] |>
    setkey(plant_id_eia)
  
# Change in load from solar panels --------------------------------------------
  panel_load_dt = 
    merge(
      nerc_state_xwalk, 
      state_solar_profile_dt, 
      by = 'state_fips'
    ) |>
    merge(
      panel_output_dt,
      by = 'nerc_adj'
    ) %>% 
    .[,.(
      nerc_adj, datetime_utc = utc_time, 
      panel_mwh = elec_profile*avg_output_mwh
    )] |>
    dcast(
      datetime_utc ~ nerc_adj,
      value.var = 'panel_mwh'
    )
  # Merging with actual load dt
  nerc_panel_load_dt = 
    merge(
      nerc_load_dt,
      panel_load_dt, 
      by = 'datetime_utc'
    ) %>% .[,':='(
      excess_load_cal_adj = excess_load_cal - cal, 
      excess_load_sq_cal_adj = (excess_load_cal - cal)^2,
      excess_load_mro_adj = excess_load_mro - mro, 
      excess_load_sq_mro_adj = (excess_load_mro - mro)^2,
      excess_load_npcc_adj = excess_load_npcc - npcc, 
      excess_load_sq_npcc_adj = (excess_load_npcc - npcc)^2,
      excess_load_rfc_adj = excess_load_rfc - rfc, 
      excess_load_sq_rfc_adj = (excess_load_rfc - rfc)^2,      
      excess_load_serc_adj = excess_load_serc - serc, 
      excess_load_sq_serc_adj = (excess_load_serc - serc)^2,
      excess_load_tre_adj = excess_load_tre - tre, 
      excess_load_sq_tre_adj = (excess_load_tre - tre)^2,
      excess_load_wecc_adj = excess_load_wecc - wecc, 
      excess_load_sq_wecc_adj = (excess_load_wecc - wecc)^2
    )]
  
# Calculating production ------------------------------------------------------
  nerc_gen_dt =
    CJ(
      plant_id_eia = model_dt$plant_id_eia,
      datetime_utc = nerc_panel_load_dt$datetime_utc
    ) |>
    merge(
      model_dt,
      by = 'plant_id_eia'
    ) |>
    merge(
      nerc_panel_load_dt, 
      by = 'datetime_utc'
    ) %>% 
    .[,':='(
        fit_net_generation_mwh_raw = intercept + 
          coef_excess_load_cal*excess_load_cal +
          coef_excess_load_mro*excess_load_mro +
          coef_excess_load_npcc*excess_load_npcc +
          coef_excess_load_rfc*excess_load_rfc +
          coef_excess_load_serc*excess_load_serc +
          coef_excess_load_tre*excess_load_tre +
          coef_excess_load_wecc*excess_load_wecc +
          coef_excess_load_sq_cal *excess_load_sq_cal +
          coef_excess_load_sq_mro *excess_load_sq_mro +
          coef_excess_load_sq_npcc*excess_load_sq_npcc +
          coef_excess_load_sq_rfc *excess_load_sq_rfc +
          coef_excess_load_sq_serc*excess_load_sq_serc +
          coef_excess_load_sq_tre *excess_load_sq_tre +
          coef_excess_load_sq_wecc*excess_load_sq_wecc,
        fit_net_generation_mwh_raw_cal = intercept + 
          coef_excess_load_cal* excess_load_cal_adj +
          coef_excess_load_mro* excess_load_mro +
          coef_excess_load_npcc*excess_load_npcc +
          coef_excess_load_rfc* excess_load_rfc +
          coef_excess_load_serc*excess_load_serc +
          coef_excess_load_tre* excess_load_tre +
          coef_excess_load_wecc*excess_load_wecc +
          coef_excess_load_sq_cal *excess_load_sq_cal_adj +
          coef_excess_load_sq_mro *excess_load_sq_mro +
          coef_excess_load_sq_npcc*excess_load_sq_npcc +
          coef_excess_load_sq_rfc *excess_load_sq_rfc +
          coef_excess_load_sq_serc*excess_load_sq_serc +
          coef_excess_load_sq_tre *excess_load_sq_tre +
          coef_excess_load_sq_wecc*excess_load_sq_wecc,
        fit_net_generation_mwh_raw_mro = intercept + 
          coef_excess_load_cal* excess_load_cal +
          coef_excess_load_mro* excess_load_mro_adj +
          coef_excess_load_npcc*excess_load_npcc +
          coef_excess_load_rfc* excess_load_rfc +
          coef_excess_load_serc*excess_load_serc +
          coef_excess_load_tre* excess_load_tre +
          coef_excess_load_wecc*excess_load_wecc +
          coef_excess_load_sq_cal *excess_load_sq_cal +
          coef_excess_load_sq_mro *excess_load_sq_mro_adj +
          coef_excess_load_sq_npcc*excess_load_sq_npcc +
          coef_excess_load_sq_rfc *excess_load_sq_rfc +
          coef_excess_load_sq_serc*excess_load_sq_serc +
          coef_excess_load_sq_tre *excess_load_sq_tre +
          coef_excess_load_sq_wecc*excess_load_sq_wecc,
        fit_net_generation_mwh_raw_npcc = intercept + 
          coef_excess_load_cal* excess_load_cal +
          coef_excess_load_mro* excess_load_mro +
          coef_excess_load_npcc*excess_load_npcc_adj +
          coef_excess_load_rfc* excess_load_rfc +
          coef_excess_load_serc*excess_load_serc +
          coef_excess_load_tre* excess_load_tre +
          coef_excess_load_wecc*excess_load_wecc +
          coef_excess_load_sq_cal *excess_load_sq_cal +
          coef_excess_load_sq_mro *excess_load_sq_mro +
          coef_excess_load_sq_npcc*excess_load_sq_npcc_adj +
          coef_excess_load_sq_rfc *excess_load_sq_rfc +
          coef_excess_load_sq_serc*excess_load_sq_serc +
          coef_excess_load_sq_tre *excess_load_sq_tre +
          coef_excess_load_sq_wecc*excess_load_sq_wecc,
        fit_net_generation_mwh_raw_rfc = intercept + 
          coef_excess_load_cal* excess_load_cal +
          coef_excess_load_mro* excess_load_mro +
          coef_excess_load_npcc*excess_load_npcc +
          coef_excess_load_rfc* excess_load_rfc_adj +
          coef_excess_load_serc*excess_load_serc +
          coef_excess_load_tre* excess_load_tre +
          coef_excess_load_wecc*excess_load_wecc +
          coef_excess_load_sq_cal *excess_load_sq_cal +
          coef_excess_load_sq_mro *excess_load_sq_mro +
          coef_excess_load_sq_npcc*excess_load_sq_npcc +
          coef_excess_load_sq_rfc *excess_load_sq_rfc_adj +
          coef_excess_load_sq_serc*excess_load_sq_serc +
          coef_excess_load_sq_tre *excess_load_sq_tre +
          coef_excess_load_sq_wecc*excess_load_sq_wecc,
        fit_net_generation_mwh_raw_serc = intercept + 
          coef_excess_load_cal* excess_load_cal +
          coef_excess_load_mro* excess_load_mro +
          coef_excess_load_npcc*excess_load_npcc +
          coef_excess_load_rfc* excess_load_rfc +
          coef_excess_load_serc*excess_load_serc_adj +
          coef_excess_load_tre* excess_load_tre +
          coef_excess_load_wecc*excess_load_wecc +
          coef_excess_load_sq_cal *excess_load_sq_cal +
          coef_excess_load_sq_mro *excess_load_sq_mro +
          coef_excess_load_sq_npcc*excess_load_sq_npcc +
          coef_excess_load_sq_rfc *excess_load_sq_rfc +
          coef_excess_load_sq_serc*excess_load_sq_serc_adj +
          coef_excess_load_sq_tre *excess_load_sq_tre +
          coef_excess_load_sq_wecc*excess_load_sq_wecc,
        fit_net_generation_mwh_raw_tre = intercept + 
          coef_excess_load_cal* excess_load_cal +
          coef_excess_load_mro* excess_load_mro +
          coef_excess_load_npcc*excess_load_npcc +
          coef_excess_load_rfc* excess_load_rfc +
          coef_excess_load_serc*excess_load_serc +
          coef_excess_load_tre* excess_load_tre_adj +
          coef_excess_load_wecc*excess_load_wecc +
          coef_excess_load_sq_cal *excess_load_sq_cal +
          coef_excess_load_sq_mro *excess_load_sq_mro +
          coef_excess_load_sq_npcc*excess_load_sq_npcc +
          coef_excess_load_sq_rfc *excess_load_sq_rfc +
          coef_excess_load_sq_serc*excess_load_sq_serc +
          coef_excess_load_sq_tre *excess_load_sq_tre_adj +
          coef_excess_load_sq_wecc*excess_load_sq_wecc,
        fit_net_generation_mwh_raw_wecc = intercept + 
          coef_excess_load_cal* excess_load_cal +
          coef_excess_load_mro* excess_load_mro +
          coef_excess_load_npcc*excess_load_npcc +
          coef_excess_load_rfc* excess_load_rfc +
          coef_excess_load_serc*excess_load_serc +
          coef_excess_load_tre* excess_load_tre +
          coef_excess_load_wecc*excess_load_wecc_adj +
          coef_excess_load_sq_cal *excess_load_sq_cal +
          coef_excess_load_sq_mro *excess_load_sq_mro +
          coef_excess_load_sq_npcc*excess_load_sq_npcc +
          coef_excess_load_sq_rfc *excess_load_sq_rfc +
          coef_excess_load_sq_serc*excess_load_sq_serc +
          coef_excess_load_sq_tre *excess_load_sq_tre +
          coef_excess_load_sq_wecc*excess_load_sq_wecc_adj
    )]
  # Drawing epsilons 
  N = nrow(nerc_panel_load_dt)
  set.seed(1234)
  epsilon_dt = 
    map_dfr(
      model_dt$plant_id_eia,
      \(id){
        data.table(
          plant_id_eia = as.character(id),
          datetime_utc = nerc_panel_load_dt$datetime_utc,
          epsilon = rnorm(
            N, mean = 0, 
            sd = exp(model_dt[plant_id_eia == as.character(id)]$log_scale)
          )
        )
      }
    ) |> setkey(plant_id_eia, datetime_utc)
  nerc_gen_dt = 
    merge(
      nerc_gen_dt[,-c('epsilon','namepcap','stack_height')],
      epsilon_dt, 
      by = c('plant_id_eia','datetime_utc')
    ) |>
    merge(
      oge_plant_info_dt, 
      by = 'plant_id_eia'
    )
  # Censoring the fitted values
  nerc_gen_dt[,':='(
    fit_net_generation_mwh = fcase(
      fit_net_generation_mwh_raw + epsilon >= 0 & 
        fit_net_generation_mwh_raw + epsilon <= namepcap, 
      fit_net_generation_mwh_raw + epsilon,
      fit_net_generation_mwh_raw + epsilon < 0, 0,
      fit_net_generation_mwh_raw + epsilon > namepcap, namepcap
    ),
    fit_net_generation_mwh_cal = fcase(
      fit_net_generation_mwh_raw_cal + epsilon >= 0 
      & fit_net_generation_mwh_raw_cal + epsilon <= namepcap, 
      fit_net_generation_mwh_raw_cal + epsilon,
      fit_net_generation_mwh_raw_cal + epsilon < 0, 0,
      fit_net_generation_mwh_raw_cal + epsilon > namepcap, namepcap
    ),
    fit_net_generation_mwh_mro = fcase(
      fit_net_generation_mwh_raw_mro + epsilon >= 0 
      & fit_net_generation_mwh_raw_mro + epsilon <= namepcap, 
      fit_net_generation_mwh_raw_mro + epsilon,
      fit_net_generation_mwh_raw_mro + epsilon < 0, 0,
      fit_net_generation_mwh_raw_mro + epsilon > namepcap, namepcap
    ),
    fit_net_generation_mwh_npcc = fcase(
      fit_net_generation_mwh_raw_npcc + epsilon >= 0 
      & fit_net_generation_mwh_raw_npcc + epsilon <= namepcap, 
      fit_net_generation_mwh_raw_npcc + epsilon,
      fit_net_generation_mwh_raw_npcc + epsilon < 0, 0,
      fit_net_generation_mwh_raw_npcc + epsilon > namepcap, namepcap
    ),
    fit_net_generation_mwh_rfc = fcase(
      fit_net_generation_mwh_raw_rfc + epsilon >= 0 
      & fit_net_generation_mwh_raw_rfc + epsilon <= namepcap, 
      fit_net_generation_mwh_raw_rfc + epsilon,
      fit_net_generation_mwh_raw_rfc + epsilon < 0, 0,
      fit_net_generation_mwh_raw_rfc + epsilon > namepcap, namepcap
    ),
    fit_net_generation_mwh_serc = fcase(
      fit_net_generation_mwh_raw_serc + epsilon >= 0 
      & fit_net_generation_mwh_raw_serc + epsilon <= namepcap, 
      fit_net_generation_mwh_raw_serc + epsilon,
      fit_net_generation_mwh_raw_serc + epsilon < 0, 0,
      fit_net_generation_mwh_raw_serc + epsilon > namepcap, namepcap
    ),
    fit_net_generation_mwh_tre = fcase(
      fit_net_generation_mwh_raw_tre + epsilon >= 0 
      & fit_net_generation_mwh_raw_tre + epsilon <= namepcap, 
      fit_net_generation_mwh_raw_tre + epsilon,
      fit_net_generation_mwh_raw_tre + epsilon < 0, 0,
      fit_net_generation_mwh_raw_tre + epsilon > namepcap, namepcap
    ),
    fit_net_generation_mwh_wecc = fcase(
      fit_net_generation_mwh_raw_wecc + epsilon >= 0 
      & fit_net_generation_mwh_raw_wecc + epsilon <= namepcap, 
      fit_net_generation_mwh_raw_wecc + epsilon,
      fit_net_generation_mwh_raw_wecc + epsilon < 0, 0,
      fit_net_generation_mwh_raw_wecc + epsilon > namepcap, namepcap
    )
  )]
  # Enforcing that production cannot increase 
  nerc_gen_dt[,':='(
      fit_net_generation_mwh_cal = ifelse(
      fit_net_generation_mwh_cal > fit_net_generation_mwh, fit_net_generation_mwh, 
      fit_net_generation_mwh_cal
    ),
      fit_net_generation_mwh_mro = ifelse(
      fit_net_generation_mwh_mro > fit_net_generation_mwh, fit_net_generation_mwh, 
      fit_net_generation_mwh_mro
    ),
      fit_net_generation_mwh_npcc = ifelse(
      fit_net_generation_mwh_npcc > fit_net_generation_mwh, fit_net_generation_mwh, 
      fit_net_generation_mwh_npcc
    ),
      fit_net_generation_mwh_rfc = ifelse(
      fit_net_generation_mwh_rfc > fit_net_generation_mwh, fit_net_generation_mwh, 
      fit_net_generation_mwh_rfc
    ),
      fit_net_generation_mwh_serc = ifelse(
      fit_net_generation_mwh_serc > fit_net_generation_mwh, fit_net_generation_mwh, 
      fit_net_generation_mwh_serc
    ),
      fit_net_generation_mwh_tre = ifelse(
      fit_net_generation_mwh_tre > fit_net_generation_mwh, fit_net_generation_mwh, 
      fit_net_generation_mwh_tre
    ),
      fit_net_generation_mwh_wecc = ifelse(
      fit_net_generation_mwh_wecc > fit_net_generation_mwh, fit_net_generation_mwh, 
      fit_net_generation_mwh_wecc
    )
  )]
  
  # Merging with damages
  tot_md_dt = 
    fread(
      here('Data/electricity-generation/output/damage-per-mwh-oge.csv')
    )[,.(
      plant_id_eia = as.character(plant_id_eia), 
      net_generation_mwh_split, 
      md_per_mwh_low, 
      md_per_mwh_high
    )] |> setkey(plant_id_eia)
  nerc_gen_dt = 
    merge(
      nerc_gen_dt[,-c('net_generation_mwh_split','md_per_mhw_low','md_per_mwh_high')],
      tot_md_dt,
      by = 'plant_id_eia'
    )
  # Caclulating damages
  nerc_gen_dt[,
    damages_baseline := 
      ifelse(
        fit_net_generation_mwh >= net_generation_mwh_split,
        (fit_net_generation_mwh-net_generation_mwh_split)*md_per_mwh_high
        + net_generation_mwh_split*md_per_mwh_low,
        fit_net_generation_mwh*md_per_mwh_low
      )
    ][,':='(
      damages_cal = 
        ifelse(
          fit_net_generation_mwh_cal >= net_generation_mwh_split,
          (fit_net_generation_mwh_cal-net_generation_mwh_split)*md_per_mwh_high
          + net_generation_mwh_split*md_per_mwh_low,
          fit_net_generation_mwh_cal*md_per_mwh_low 
        ),
      damages_mro = 
        ifelse(
          fit_net_generation_mwh_mro >= net_generation_mwh_split,
          (fit_net_generation_mwh_mro-net_generation_mwh_split)*md_per_mwh_high
          + net_generation_mwh_split*md_per_mwh_low,
          fit_net_generation_mwh_mro*md_per_mwh_low 
        ),
      damages_npcc = 
        ifelse(
          fit_net_generation_mwh_npcc >= net_generation_mwh_split,
          (fit_net_generation_mwh_npcc-net_generation_mwh_split)*md_per_mwh_high
          + net_generation_mwh_split*md_per_mwh_low,
          fit_net_generation_mwh_npcc*md_per_mwh_low 
        ),
      damages_rfc = 
        ifelse(
          fit_net_generation_mwh_rfc >= net_generation_mwh_split,
          (fit_net_generation_mwh_rfc-net_generation_mwh_split)*md_per_mwh_high
          + net_generation_mwh_split*md_per_mwh_low,
          fit_net_generation_mwh_rfc*md_per_mwh_low 
        ),
      damages_serc = 
        ifelse(
          fit_net_generation_mwh_serc >= net_generation_mwh_split,
          (fit_net_generation_mwh_serc-net_generation_mwh_split)*md_per_mwh_high
          + net_generation_mwh_split*md_per_mwh_low,
          fit_net_generation_mwh_serc*md_per_mwh_low 
        ),
      damages_tre = 
        ifelse(
          fit_net_generation_mwh_tre >= net_generation_mwh_split,
          (fit_net_generation_mwh_tre-net_generation_mwh_split)*md_per_mwh_high
          + net_generation_mwh_split*md_per_mwh_low,
          fit_net_generation_mwh_tre*md_per_mwh_low 
        ),
      damages_wecc = 
        ifelse(
          fit_net_generation_mwh_wecc >= net_generation_mwh_split,
          (fit_net_generation_mwh_wecc-net_generation_mwh_split)*md_per_mwh_high
          + net_generation_mwh_split*md_per_mwh_low,
          fit_net_generation_mwh_wecc*md_per_mwh_low 
        )
  )]

  # Aggregating over all power plants 
  nerc_gen_dt[,.(
    cal = sum(damages_cal- damages_baseline),
    mro = sum(damages_mro- damages_baseline),
    npcc = sum(damages_npcc- damages_baseline),
    rfc = sum(damages_rfc- damages_baseline),
    serc = sum(damages_serc- damages_baseline),
    tre = sum(damages_tre- damages_baseline),
    wecc = sum(damages_wecc- damages_baseline)
  ), by = .(shaped = as.numeric(plant_id_eia) > 900000)] |>
    melt(
      measure.vars = 2:8,
      variable.name = 'nerc_adj',
      value.name = 'damages'
    ) |>
    dcast(
      nerc_adj ~ shaped, value.var = 'damages'
    ) %>% .[,.(
      nerc_adj, 
      damages = `TRUE` + `FALSE`,
      Normal = `FALSE`,
      Fleets = `TRUE`
    )]
  # Checking the predicted change in production 
  nerc_gen_dt[,.(
    cal  = sum(fit_net_generation_mwh_cal - fit_net_generation_mwh),
    mro  = sum(fit_net_generation_mwh_mro - fit_net_generation_mwh),
    npcc = sum(fit_net_generation_mwh_npcc - fit_net_generation_mwh),
    rfc  = sum(fit_net_generation_mwh_rfc - fit_net_generation_mwh),
    serc = sum(fit_net_generation_mwh_serc - fit_net_generation_mwh),
    tre  = sum(fit_net_generation_mwh_tre - fit_net_generation_mwh),
    wecc = sum(fit_net_generation_mwh_wecc - fit_net_generation_mwh)
  ), by = .(shaped = as.numeric(plant_id_eia) > 900000)]|>
    melt(
      measure.vars = 2:8,
      variable.name = 'nerc_adj',
      value.name = 'change_in_production'
    ) |>
  merge(
    panel_output_dt,
    by = 'nerc_adj'
  )  
# Issue: Don't match with Mark, Why? ---------------------------------------
  p_load(readxl, janitor, ggplot2)
  mark_results = read_xlsx(here('StructuralCode/PPProdYearBaseline2.xlsx')) |>
    clean_names() |>
    data.table() %>% 
    .[-1,.(
      baseline_production_by_pp = 16*as.numeric(baseline_production_by_pp), 
      add_one_panel_in_ohio = 16*as.numeric(add_one_panel_in_ohio),
      mark_baseline_damages = 16*as.numeric(damages_3)*1000, 
      mark_ohio_damages = 16*as.numeric(damages_19)*1000
    )]
  pp_dt = nerc_gen_dt[,.(
    baseline = sum(fit_net_generation_mwh),
    prod_ohio = sum(fit_net_generation_mwh_rfc),
    damages_baseline = sum(damages_baseline),
    damages_ohio = sum(damages_rfc)),
    keyby = .(plant_id_eia=as.numeric(plant_id_eia))
  ]
  pp_dt = 
    cbind(
      pp_dt, 
      mark_results
    )
  pp_dt[plant_id_eia !=2410,.(
    tot_change_prod_emmett =  sum(baseline) - sum(prod_ohio),
    tot_change_prod_mark =  sum(baseline_production_by_pp) - sum(add_one_panel_in_ohio),
    baseline_damages= sum(damages_baseline),
    ohio_damages =  sum(damages_ohio),
    tot_change_damages_emmett =  sum(damages_baseline) - sum(damages_ohio)
  )]
    ggplot(pp_dt, aes(
      x = prod_ohio-baseline, 
      y = add_one_panel_in_ohio-baseline_production_by_pp)
    ) +
      geom_point() + 
      geom_abline(slope = 1 , intercept = 0, linetype = 'dashed') + 
      #ylim(-0.12,0.015)  +
      #xlim(-0.12,0.015) +
      labs(x = 'Emmett Change in Production', y = 'Mark Change in Production') + 
      theme_minimal()
    
    
    ggplot(pp_dt[plant_id_eia != 2410], aes(x =damages_ohio-damages_baseline , y = mark_ohio_damages-mark_baseline_damages)) +
      geom_point() + 
      geom_abline(slope = 1 , intercept = 0)
  
  
  
  